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This was supposed to be the Edition 1 of the #3 report presented to Professor Katsuhisa Furuta 
by Kittisak Tiyapan. Then I quitted my Ph.D. degree study. So the technical report became a 
technical record. The number still remained the same. Upto this moment there had already been 
three technical reports, namely Technical Report #1, #2, and #6. From now on I will call the 
rest by the name Technical Record. This is also good in a way that it will not be mistaken for a 
typical technical report published by an instution. I have no intention to claim this as a published 
item when it is not. 

All simulation for this report was done by MATLAB Version 5.1.0.421 on a Solaris 2 machine 
named Hayate. 


Tunnel diode circuit 


Consider a Tunnel Diode circuit of Figure 1 [Kha96] 


L 


Figure 1: Tunnel Diode Circuit. 


which in state space model is represented by 


[—A(21) + X2| 


1 = 


L2 = [—21 — Rro+ul. 


BIE Q)]R 
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(1) 
(2) 


Equation 1 and 2 are based on the assumption that both the capacitor C’ and the inductor L are 


linear and time-invariant. Their models are 


duc dir, 
= = |,— 
IE gS 
The tunnel diode is characterized by 
tr = h(vr) 


Then by Kirchhoff’s current law, 
io tig—-ip=0 3S ig=—h(a1)4+ 22. 
And by Kirchhoff’s voltage law, 
ve —-E4+Rip+tv,=0 => vp, = 4, — Rrot+u. 
And where t1=U0, ta=i,, u=E. 


The equilibrium point correspond to the roots of 


E 


1 
h (x1) = R = R 


Let the areu=3V, R=5kOQ, C=6 pF and L =7 wH. Then the state-space model with 


time measured in nanoseconds is (the current x2 and h(21) are in mA) 


[—h(x1) + £9] 


t1 = 


me OlrR 


2 = =[-2% —52_.+3]). 
Let the characteristic of this diode be represented by 
h(xi) = 20x, — 100x? + 200z} + 200x} + 80x? 
which is graphically shown as Figure 2 (a). Then let the characteristic be 


h(x1) = 182, — 104x7 + 230x} — 230x} + 842} 


(6) 


which is Figure 2 (b). Figure 2 shows that the characteristic of the Tunnel Diode for Equation 5 


and Equation 6 look similar graphically. 
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Tunnel giode chayacterstic, 5 
x10° Tunnel Diode characteristic n= 18 Vv, — 104 ve + 230 ve - 230 vi, + 84 ve 
T T T 


Voltage 


Figure 2: Tunnel Diode characteristic 


Some knowledge about this system 


To find equilibrium points obtain Equation 7 and 8 from Equation 3 and 4 respectively. 


1 
ae [—h(x1) + £9] (7) 
1 
0 = 7 [#1 — 522 + 3). (8) 
Solving Equation 7, 8, and 5 leads to the equilibrium points 
—3 — x 
(eZ): 9) 
where “a, is the x, at a equilibrium point and can be obtained by solving the equation 
—3 + 1012, — 5002? + 1000 x? + 1000 x} + 4002? = 0. (10) 
or the equation 
eee 
182, — 104a? + 2302? — 23024 + 84a? — preci (11) 


for the case when the characteristic of the diode is Equation 5 and 6 respectively. The solution to 


Equation 10 is 
¢, = —1.4541.221, 0.0355, 0.184+0.1572, 


while that for Equation 11 is 
x, © 0.042257, 0.314562, 0.603715 + 0.3054667, 1.173846. 


In the first case the only equilibrium point is the point where 2 is a real root, ie (0.0355, —0.6071), 
while in the second case there are three equilibrium points, ie (0.042257, —0.6085), (0.314562, —0.6629), 
and (1.1738, —0.8348). 

From Equation 3 to 5 the state equation of the system can be written as 


= file = ; (—2021 + 100x? — 200x3 — 200x¢ — 802° + 2x2) 


fg = folz) = 7[-21 — 522 +3] a) 
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Find the Jacobian matrix from 


Ofi(z) Afi (x) 
7 = | fate) Oats | (13) 
we 0x1 Ox2 
and 
Ofi(z) _ —20a, + 100x? — 200x3 — 200x} — 8027+ 22 Of; (x) _1 Ofx(t)_ 1 Ofi(z)_ 5 
Ox, 6 , Ox 6’ Ox, T Ox iC 


which leads to 


—0.1429 —0.7143 


7 —2.2821 0.1667 | 


From which the eigenvalues can be obtained to be at —2.2668 and —0.7296. Both are negative, 
therefore are stable nodes. 
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Figure 3 shows the phase plane of the diode of Equation 3 to 5. The initial conditions used in 
simulating Figure 3 was {{—5,5} x {-5 <Z < 5}} U{{-5 <Z< 5} x {-5,5}}. 


Phase plane of a Tunnel Diode 


x2 
{= | 


1 1 i it i 
5 -4 -3 -2 -1 0 1 2 3 4 5 


Figure 3: Phase plane, Tunnel Diode. 


With a sign function in feedback 
Figure 4 shows the phase plane of the diode when u = = + sgn(x, + X2). 
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Phase plane of a Tunnel Diode 


} 
-1 
-2 
-3 
-4 
5 1 1 1 = i 1 i 
-3 -2 -1 0) 1 2 3 
x1 


Figure 4: Phase plane with sign state feed back, Tunnel Diode. 


x2 
oO 
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With a sign function and state variables in the feedback 


Figure 4 shows the phase plane of the diode when u = 2 + sgn(x1 + x2) — (x1 + Xa). 


Phase plane of a Tunnel Diode 


x2 
lo} 


Figure 5: Phase plane with sign and state feed back, Tunnel Diode. 


Backstepping control design 


Consider the case where the characteristic of the diode is not known, then h(x) is uncertain, but 
is in the form of a fifth order polynomial passing through the origin (ie in the form of Equation 
5 and 6 with coefficients uncertain but with bounds.) Then Equation 3 and 4 do not satisfy the 
matching condition since there is an uncertain term in Equation 3 where there is no control input. 
The resistance R can vary. The capacitance C’ and the inductance L are assumed to be exactly 


known. In other words, the system can be written as 


Ly = 0,24 Ox” O32°3 0425 Osx° + 8X9 
tg = —rx,+ 6x2 +714, 
where |01| <a, |@2| < b,|03| < c, |O4| < d,|05| <e,06 <p 
and a, b, C, d,€,D, r, 8, O6 S Re. 


Notice that when the characteristic of the tunnel diode is described by 


2 3 4 5 
h(ai) = K1%1 + Kot] + Kx} + Kar] + 52}, 


Let the parameters be as described in page 2 and h(21) as in Equation 5, the result is shown in 
Figure 6. Parameters used for simulating Figure 6 may be summarized to comply with Equation 
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a 3 Tunnel diode 
h(x,)= 20 x,- 100 x}+ 200 x1+ 200 x/+ 80 x?, u= 3V, R=5 kQ, C= 6pF, L =7uH 


1° ; f T T T 


Figure 6: Tunnel Diode, h(x,) = 20x,— 100x7+ 200r$+ 200x{+ 80x}?, u = 3V, R = 5kO, C = 6pF, 
L=7pH 
14, 15 and as 

tq =, 20K = —100,. he = 200, 64 = 200, np = 80, BH 38V, R= bkO, C]btpk, L=T pd, 


time is measured in nanoseconds and current in mA. The time of simulation for Figure 6 is 5 seconds. 
The initial conditions are Initial Condition Set 1. 


Initial Condition Set 1 The set of points (x1(0),x2(0)) created by 


n = QO; 
for 1 = -1 to 1 steps 0.3 
n=n+1, x1(0) =i and x2(0) = i are the n_th initial conditions 
n=n+1, x1(0) = i and x2(0) = -i are the n_th initial conditions 
endfor 


Then let the parameters be the same but let h(2,) be as in Equation 6 instead, the result is 
shown in Figure 7. Parameters used for simulating Figure 7 may be summarized to comply with 
Equation 14, 15 and as 


Ky =,18 Ko = —104, K3 = 230, Ka = —230, K5 = 84, E = 3V, = 5.6: C = 6pF, = Cid. 


The initial conditions were Initial Condition Set 1. The time of simulation is 5 seconds for half of 
the initial conditions set (the half with (7,7) according to Initial Condition Set 1) and 10 seconds 
for the rest (ie the half with (i, —7).) 

Suppose that the uncertainties within this system are described as follows, 


—-5 < ki < 30 ; —200 < ke < O 
10 < «3 < 500 ; —500 < Ka <_ 500 
S10" rope “LOO ’ LQ: dS TORO 
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h(x,)= 18 x,— 104 x24 230 x9 230 eeean u= 8V, R=5 kQ, C= 6pF, L =7uH 

19 - 1 1 
0.8F 4 
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Figure 7: Tunnel Diode, h(x,) = 182,— 104x%?+ 230x73— 230xr1+ 842}, u = 3V, R = 5kQ, C = 6pF, 
L= pH 


which could be rewritten (though with somewhat more stringent a requirement) as 


al < 8 , [Os] < 20 , [Ol < M=1.4286 


The following theorem is actually Lemma 13.3 in Khalil [Kha96]. 


Theorem 1 Consider the system 


a = Fm) + 9m + n(n, €) (16) 
€ = faln,€) + ga(n, €)u + de(n, €), (17) 
defined on a domain D C R"*" that contains the origin. Let the uncertainty satisfies V(n, €) € D, 
onl, Olle < erlinll. (18) 
e(7,€)| << aa [InIlp + asl é| (19) 
Let o(n) be a stabilizing state feedback control for Equation 16 that satisfies 
1 
Val. €) = V(n) + 5[E - on)? (20) 
and V(n) be a Lyapunov function that satisfies Vn € D 
O 
jail <ul — [5° <a (21) 
Then the state feedback control 
1 |0¢ 


u=— 


OV 
- Sor +6) - Fa fe HE- 9), an (22) 
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with k sufficiently large, stabilizes the origin of Equation 16 and 17. Moreover, if all the assumptions 
hold globally and V(n) is radially unbounded, the origin will be globally asymptotically stable. a 


(see Khalil [Kha96] Lemma 13.3 for the proof) 
From Theorem 1 Equation 14 and 15 is of the form of Equation 16 and 17 respectively, where 


1) = Ty ) c = wr 
f(n) = 0 , gm) = 8 
On = O12, + Oox7 + Ogz$ + Gazi + Osx? , faln,€) = rx 
ga(n, §) ames ; dg = x2 


The uncertainty terms are 


64 = On = A424 Ox" 30° Oar O5xr°, and 9 = de = AgX2. 


For |a| <¢ and |r2| <o 


|51] < 56 +356? + 85¢% + 85¢4 + 206° 
< (5+ 35¢ + 85<? + 855% + 20¢4) |x| 
and |5p| < 1.42860 < 1.4286 |zr9). 
Consider Equation 14 as a system and view x2 as a control input (notice that s = é pS 7) 
r2> o1(#1) = —kyx, 
1 
and Vil) = 51 
Then 2 
Vi = 8x10(x1) + 0127 + Orr} + O32} + O42} + Osx} 


< — (sk; — (5 + 35¢ + 85c? + 85<% + 20¢*)) af. 


Choose , — 1 +41 +25 + O35" + Bas” + O55" _ 1+ 5 + 35s + 855° + 85s" + 20c" 
1 ; : : 
Then from Equation 20 and 22 


x 1 1 
— 7|-n 2 = (45) —k (ea +hn)| 


= —1.1669k x2 — 2.16652, — 0.1429k (aq + ky21), and 
Va 


; [xt + (r+ kye1)"| 


The derivative of the Lyapunov function is 


. 1 
V, = 3 [2141 + 2 (kyary + x2) (kiay + £2) 
Restrict the analysis to the set 
Ope {ce R | Va(x) Kel. 


Notice that 9 
es 1 c | (« _ (1+ max |6;| + max |92|¢ + max |43|¢? + max |44|¢° toe Pel | 
= 9 l l by 


S 


Substituting 
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= 042, + Oor? + O32? + O4r + Osx? — ky sx 
fg = —2rax, — sx, — kysxq + O6%2 —k (kya, 4+ 22) 
gives 
V, = ; [22° (-skir1 + O01 + Oor} + O3r3 + O4xt 4 O52) 
+2 (kx, + £2) (-2ra1 — sx, + ky (-skixi + O12 + Ooay + 0323 + O4x} 4 O53) 
— sky tq + 602 — k (kyx1 + £2))| 
= 22 + (kyr, +22) (-2r —s- kis) xy + (kyr, + 22) Gi + Oox + O30? + O42 4 O54) ee Ry 
+ (kx, + rq) (—sk, + 06) 22 — k (kya, + at)" 
< (kya, + 22) (-2r —s- kis) 21+ |kyxy + £9! (1 t Aa¢ + 0367 + O45? 4 5c") ¢ (kix1 + £2 — Le) 
+ (kx, + £2) (—sk, + 06) vo — k (kya, + at)" 
< (kya, + 22) (-2r —s- kis) x1 + |kyx, + £9! (-1 —(1+¢) (1 + Oo¢ + 0357 + O45? 4 5c") c+ 06 
_ (k - (1 + B25 + 0367 + O45? 4 65<*) <) (kyx1 + a2)” 
Choose 


2 
k > (max |@1| + max |O2|¢ + max |03|¢? + max |@4|c® + max \05 |<") ot+ (—2r —s— kis) 


2 
+ [max [06] -1 —(1 +9) (max |@1| + max |05|¢ + max |03|¢? + max |94|c? + max \4s|s*) | 


would ensure that the origin is exponentially stable (see Khalil [Kha96] Corollary 3.4), with Q, 
contained in the region of attraction. 
Then with the values of parameter above and 


¢=0.1, 0=0.2 


the gain k becomes 
ky = 44-62 SF 1,105 10". 


To show graphically the result, consider the system where 
Ky = 22, Kg = —180, K3 = 440, K_499 =, K5 = 90, R= LV, R= 9kQ, C= 6 pF, L= 7 LA, SoS 0.1, Oo 
The initial conditions used are Initial Condition Set 2 and the simulation time is 5 seconds. 


Initial Condition Set 2 The set of points (x1(0),22(0)) created by 


n = 0; 
for 1 = -0.5 to 0.5 steps 0.11 
n=n+1, x1(0) =i and x2(0) = i are the n_th initial conditions 
n=n+1, x1(0) = i and x2(0) = -i are the n_th initial conditions 
endfor 


The result is shown in Figure 8 Apply a control input as shown in Equation 23, k,,k as has been 
calculated above. The simulation result is shown in Figure 9 (Initial Condition Set 2 and simulation 
time 5 seconds.) 
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Tungel diode 
h(x,)= 22 x,— 180 x24 440 x3_ 490 x,+ 90 x), u= 1V, R=9 kQ, C= 6pF, L =71H 
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Figure 8: Tunnel Diode, h(x1) = 22x%1— 18027+ 440r3— 490x{+ 9027, u=1V, R= 9kO, C = 6pF, 
C= Tin 


Tynnel diode with, control 
h(x,)= 22 x,— 180 x2 + 440 x)— 490 x/+ 90 x7, u= 1V, R=9 kQ, C= 6pF, L =7uH 


25 T T T T T T 


-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 


Figure 9: Tunnel Diode with control, h(x1) = 22x,— 180x7+ 44023— 490r{+ 9022, u=1V, R= 
9kQ, C = 6pF, L = TH 
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Lorenz attractor 


13" May 1998 
Consider a Lorenz attractor 


L1 = (xq = r1) (23) 
Lo => Ly (1 + v\— r3) — % (24) 
£3 = UML — bx3. (25) 


This system has got three equilibrium points, one of which is the origin and the other two are 
symmetrically positioned with respect to each other. To find the equilibrium points, let ¢ = 0 
which leads to 


= 0 (2-21) (26) 
= (1 +X — 23) — X92 (27) 
= Wr — bx3, (28) 


from which the equilibrium points are 


(01,2203) = (0,00), (-YOVO.-VOVO.a), (YOVO. VoOVved.>) 


The equilibrium point at the origin is unstable and the other two will be stable if the condition 


A(b+1—e)+(6+1)(o+6+1)>0 (29) 


is met. 


Simulation result of the system 


Results from simulations done on this system are shown from Figure 10 to Figure 12. Figure 10 is 
the system with 
=a x=] bo: ba, 


(Equation 29 is satisfied) and the initial conditions were 


1. for (ia te dhs ele) 
2. for (b), (#1, 22,03) € {-5 <Z* <5} x {1} x {-5, 5} U {5,5} x {1} x {-5 < It <5}, 
( 


( 
( 

3. for (c), (v1, £2,%3) € {-5 < Zt <5} x {-5,5} x {1} U{—5, 5} x {-5 < Zt <5} x {1}, 

4. for (d), (w1,22,23) € {1} x {-5 < Zt <5} x {-5, 5}U {1} x {-5, 5} x {-5 < T+ <5}, and 
( 


5. for 


The CPU time in the case of Figure 10(e) was 1.07 seconds, and its file size is 98.88 KB. 
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(init.cond. (7,7,-7) to (-7,-7,7)) (init.cond. (—7,7,7) to (7,-7,-7)) 


Figure 10: Lorenz attractor (o = 3, X= 5, b= 1), (a) State variables, (b) State plane x, v.8. £3, 
(c) State plane x, v.8. £2, (d) State plane x2 v.s. £3, (e) State space 
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Figure 11 was from the system with 


L. for (eines .ag = (Led), 
2. for (b), (a1, 22,23) € {-5 < Z* < 5} x {-5, 5} x {1} U{-5,5} x {-5 < Zt < 5} x {1}, 
( 


( 

( 
3. for (c), (a1,#2, 3) € {1} x {-5 < It <5} x {-5, 5}U {1} x {-5,5} x {-5 < It <5}, 
A. for (d), (#1, 2,23) € {-5 < I+ <5} x {1} x {-5, 5}U{-5, 5} x {1} x {-5 < I+ <5}, and 
5. for ( 


The CPU time used in simulating Figure 11 (e) was 2.32 seconds, and its file size is 217.24 KB. 
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Figure 11: Lorenz attractor (o = 15, X = 25, b= 3), (a) State variables, (b) State plane x v.s. £2, 
(c) State plane x2 v.s. £3, (d) State plane x, v.s. £3, (e) State space 
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Figure 12 shows the result from a system having 


(Equation 29 is not satisfied) together with the initial conditions as shown in page 15. The CPU 
time used during the simulation for Figure 12 (e) was 2.01 seconds, and its file size 196.51 KB. 
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Figure 12: Lorenz attractor (o = 9, \ = 30, b= 1), (a) State variables (b)State plane x, v.s8. £2, 
(c) State plane x2 v.s. £3, (d) State plane x, v.s. £3, (e) State space 
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Lorenz attractor with control 


Add an input into Equation 25 so that the system is now 


Ly = C6 (xo = 2X1) (30) 
Xo = Vy (1+ A— 23) — %@ (31) 
3 = U4X_— brig + u. (32) 


To demonstrate the effect of a sign function on the response let 
u = Ksgn (x; + Xo + x3). (33) 


(Errata: the initial condition labels written on Figure 10 should be the same as that of Figure 14, 
ie. all 7’s should be replaced by 8’s) 
Let K = —1 and the input now becomes 


u = —sgn (x; + x2 + X3). (34) 
Simulated result with this input is shown in Figure 13. Here the conditions were 
G=3, = 5. b=, 


(Equation 29 is satisfied) and the states at the start of simulation were the same as those initial 
conditions shown in page 15. Figure 13 (c) has got tcpy = 19.96 seconds and the size of the file 
containing it is 1.87 MB. 
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Figure 13: (a) State variables of a Lorenz attractor (a 
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3, A = 5, b = 1), (a) State variables, 


(b) State plane x1 v.8. %, (c) State plane x2 v.s. x3, (d) State plane x, v.s. x3, (e) State space 


u = —sgn (x1 + x2 + x5) 


Technical record #3 by Kittisak Tiyapan, Tokyo Institute of Technology 2A 


Let K = —3 and obtain 
u = —3 sgn (x1 + Xo + Xs). (35) 


The conditions of the system was 


(Equation 29 is satisfied) and the initial conditions used were the same as those shown in page 15. 
The result is show in Figure 14. Figure 14 has got topy = 26.25 seconds and its file is 2.38 MB in 
size. 
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Figure 14: Lorenz attractor (o = 3, \ = 5, b = 1), tsimutation = 50 sec) (a) State variables, 
(b) State plane x, v.s. £2, (c) State plane x2 v.s. x3,(d) State plane x, v.s. x3,(e) State space. 
u = —3 sgn (x; + X2 + x3) 
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Let K = —5 and obtain 
u = —5 sgn (x; + Xo +X3). (36) 


The conditions of the system was 


(Equation 29 is satisfied) and the initial conditions used were the same to those shown in page 15. 
Figure 15 shows the result. The file size of Figure 15 (e) is 2.38 MB, and it took topy = 26.25 
seconds. 
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(e) 


Figure 15: (a) State variables of a Lorenz attractor (o = 3, X= 5, b= 1), (b) State plane x, v.s. x2 
, (c) State plane xq v.s. x3, (d) State plane x, v.s. x3, (e) State space. u = —5 sgn (x, + Xo + x3) 
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Figure 16 shows the result when the system had 


(ie. the conditions of Equation 29 was not satisfied) simulated with the same set of initial conditions 
to those shown in page 15. The input was taken as 


u = —sgn (x; + X29 + x3). 


For Figure 16 (e) tepy = 2.59 seconds, and file size = 203.62 KB. 
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Lorenz attractor with input 
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Figure 16: A Lorenz attractor (o =9, A = 30, b= 1), (a) State variables, (b) State plane x1 v.s. 
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Lorenz attractor with input 
(init.cond. (8,-8,8) to (-8,8,-8)) 
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Lorenz attractor with input 
(init.cond. (-8,8,8) to (8,-8,-8)) 


X2, (c) State plane xq v.s. £3, (d) State plane x, v.s. x3, (e), u= —sgn (x; + X2 + x3) 
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Let the input be 
u = —5 sgn (x1 + X2 + x3), 


and the system has 
=O) AS 20, “b= 1, 


so that the conditions of Equation 29 was not satisfied. Then with the set of initial conditions 
specified in page 15 simulated and obtained the result of Figure 17. Figure 17 has tepy = 2.62 
seconds, and file size = 209.80 KB. 
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Figure 17: A Lorenz attractor (o =9, A = 30, b= 1), (a) State variables, (b) State plane x1 v.s. 
X2, (c) State plane xq v.s. x3, (d) State plane x, v.s. £3, (e), w= —5sgn (x; + x2 + x3) 
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Feeding back u = k sgn(x; + xg + x3) to the system with 
k>0 


led to the response shown in Figure 18. Here the two equilibrium points away from the origin are 
not stable. Figure 18 was from a simulation with 


u = 3 sgn (x; + X2 + x3) 


as an input. The system had 
G=9.. X= 30,, .b=<1, 


so that the conditions of Equation 29 was not satisfied. Then simulated with the set of initial 
conditions the same with those given in page 15. Figure 18 (e) has tcpy = 2.64 seconds, and file 
size = 206.55 KB. 
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Figure 18: A Lorenz attractor (o =9, A = 30, b= 1), (a) State variables, (b) State plane x1 v.s. 
Xo, (c) State plane xq v.s. £3, (d) State plane x, v.s. 43, (e), u= 3 sgn (x; + Xo + x3) 
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Discussion 
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Synchronous machine 


11 May 1998 
Consider a model of a synchronous electrical machine [Coo86] used in power system engineering. 


Hd = Py, — P.sind (37) 


where 0 is the rotor angle relative to a reference frame, H is the inertia constant, P,, is the mechanical 
power supplied, and P, is the maximum electrical power which can be generated. This equation 
can be put into state-space form as 


L1 = (38) 
: u— P.sinx 

2. om (39) 
Y= P.sin x1 (40) 


where 21 = 6, 2 = 6 andu= Re. 
Figure 19 shows the state variables and the phase plane of this machine respectively when 


Bad, Pea eS. 


Phase plane 
T 


ssesaasenaee (©) 


5 6 - 
Time in seconds x1 


(a) (b) 


Figure 19: A synchronous machine (H =1, Pn =1, P. = 2), (a) State variables, and 
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Unless the mechanical power is within the bounds defined by the maximum electrical power the 
system can not operate in equlibrium. Figure’s 20 (a) and (b) shows the state variables and the 
phase plane of this machine respectively when 


Hei, Pe=3- Po: 


that-is Pp, > P.. 


Phase plane 


5 6 ae i“ 
Time in seconds x1 


Figure 20: A synchronous machine (H = 1, Py, = 3, P. = 2), (a) State variables, and (b) Phase 
plane. 
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When 
Py <0 


the machine takes power from the bus and acts as a motor instead of a generator. Figure 21 (a) 
and (b) shows the state variables and the phase plane of this machine respectively. Here 


HA Pe h.. Peso 


Phase plane 


5 6 
Time in seconds 


Figure 21: A synchronous machine (H =1, Pn = —1, P. = 2), (a) State variables, and (b) Phase 
plane. 
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If 
|Pn| > |Pe| 


the synchronous motor will not operate at equilibrium as shown in Figure 22. 
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Figure 22: A synchronous machine (H = 1, Py = —3, P. = 2), (a) State variables, and (b) Phase 
plane. 
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With a signum function feedback 


Fed the sign function of state variables back to input of the machine and obtained Figure 23. Here 
Hk. Gees. Pees. 


(a generator) 


Phase plane 
T 


5 6 - 
Time in seconds x1 


Figure 23: A synchronous machine (H = 1, P, = 1, P. = 2), (a) State variables, and (b) Phase 
plane, with a sign function feedback. 
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Figure’s 24 (a) and (b) shows the results, also with a sign function in feedback, but this time 
Pee U. 


(a motor) The system had 
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Figure 24: A synchronous motor (H =1, Pn = 
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(b), (d), (f) are Phase planes, with a sign function feedback. 


—1, P.=2), (a), (c), (e) State 


variable plots, and 
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When |P,,,| > |P.| 


Again, P,, must satisfy the condition 
\Pn| > |Pe| 


or the results similar to that shown in Figure 25 (a) and (b) will be obtained. 


Phase plane 
T 


5 
Time in seconds 


(a) (b) 


Figure 25: A synchronous motor (H =1, Py = 3, P. = 2), (a) State variable plots, and (b) Phase 
planes, with a sign function feedback. 
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Next figure, Figure 26 shows the phase plane of a condition similar to that of Figure 25 but 
where P. = —2 instead of 2. 


Phase plane 


20 T T i T 


10 15 20 


Figure 26: A synchronous motor (H = 1, Py, = 3, P. = —2), (a) State variable plots, and (b) 


Phase planes, with a sign function feedback. 
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With P,, > 0 but P. <0 


Figure 27 shows the result when conditions are similar to those of Figure 23 except that this time 
P. = —2 instead of 2. 
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Figure 27: A synchronous motor (H = 1, Pn = 1, Pe = —2), (a) State variable plots, and (b) 
Phase planes, with a sign function feedback. 
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Figure 28 (b) is similar to Figure 28 (a) except that here there is no feedback of sign function. 
By comparing the two figures one can see that about half of the equilibriums were lost when a sign 
function FB was introduced. The parameters used in simulating Figure 28 were 


Hat) Peds piso. 
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Figure 28: A a synchronous machine (H =1, Pn =1, Pe = 2). (a) with a sign function FB and, 
(b) without a sign function FB. 
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Figure 29 shows the result when 


and 
u = —Ddsgn(x, + X2) 
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Figure 29: Phase plane of a synchronous machine (H =1, Pn =1, Pe = 2). u = —5sgn(x, + x2) 
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sgn (Fin) | sgn (Fe(9)) | sgn((0)) | Figure # 

+ + Figure 30 
+ + — Figure 31 
+ — + Figure 32 
+ — = Figure 33 
— “+ + Figure 34 
— + = Figure 35 
_ -_ + Figure 36 
= _ _ Figure 37 


Table 1: Sign table for the investigation to study the effect of initial conditions on response 


Initial conditions # || |P.(0)| |5(0)| 
# 1 sgn (Fm)-1 | sgn(P.(0)) + 1.5 | sgn (6(0)) - 5 
# 2 sen (Pm)-1.5 | sgn(P.(0))-1 | sgn (6(0)) - 5 
#3 sgn (Pm) +3 | sgn(P.(0)) +3 | sgn (5(0)) 
# 4 sen(Pn)-5 | sgn(P.(0))-5 | sgn (6(0)) + F 
ae5 sgn (Pn) - 0.1 | sgn (P.(0)) -0.1 | sgn (6(0)) - § 


Table 2: The initial conditons to study the effect of initial conditions on response 


Three-dimensional model of synchronous machine 


Consider the three-dimensional model of the synchronous machine (see also |Tiy98]) described by 


Ly = (41) 
: ce ©. : 
a a We _ nts sin 21 (42) 
E 
te = a ss +B cos ce peceeaeiees (43) 
T T T 


where 


a, =6(t), «.=d(t), 23 = P(t). 


The following study follows the initial conditions and parameter set described by Table 1 and 
Table 2. The rest of the parameters are 


D 
Pry = 1, H = 0.03, = 2.2, = 2.9,93 = 17,7 = 8, Epp = 1.45, 5 = -6, 


and another initial condition still to be mentioned is 
6(0) =0. 


The following are the figures mentioned in Table 1. 
The result from Figure 30 to Figure 37 can be summarized as follows. The system was unstable 
(ie. 6 increased indefinitely) when 


1. Ph = +5, P.(0)=+5, 6(0) =+% (every figure from Figure 30 to 37), 
2. Pr =+3, P.(0)=+3, 6(0) = +% (every figure from Figure 30 to 37), 
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3. Pm =+1, P.(0)=—1.5, 6(0) =+5 (Figure 32, 33, 36, and 37), or 
4. Pp =+1.5, P(0)=—1, 46(0) =+ 5 (Figure 32, 33, 36, and 37). 
The system was stable (ie. 6 did not increase indefinitely) when 
1. Pm =+0.1, P(0)=+0.1, 6(0) =+§ (every figure from Figure 30 to 37), 
2. Pm = +1, P.(0)=1.5, 6(0) =+% (Figure 30, 31, 34, and 35), or 
3. Pm = 1.5, P.(0)=1, 46(0) = +5 (Figure 30, 31, 34, and 35). 


Notice that P. > 0 denotes electrical power produced from the machine, while P. < 0 denotes the 
opposite, ie. the electrical power injected into the machine. Therefore, to put it another words is 
that no matter whether the machine is having P,, > 0 (mechanical power input into the machine) 
or P,, < 0 (mechanical power produced by the machine), when ever P. < 0 the machine becomes 
less robust to the initial states of the state variables. Or we can say that the domain of attraction 
in the phase plane when P. < 0 is smaller when that when P, > 0. 
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‘Synchronous machine, x,(t) 
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Figure 30: 6(0) > 0, 6(0) > 0, P.(0) > 0. 
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‘Synchronous machine, x,() Synchronous machine, x,() Synchronous machine, ,() 
P,,(0)>0.P (0)>0,8(0}<0 P,,(0)>0.P,(0)>0,8(0}<0 P -(0)>0,P,(0)90,8(0)<0 
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Figure 31: 6(0) > 0, 6(0) > 0, P.(0) > 0. 
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Figure 32: 5(0) > 0, 6(0) > 0, P.(0) > 0. 
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‘Synchronous machine, x,(t) 
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Figure 33: 5(0) > 0, 6(0) > 0, P.(0) > 0. 
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Synchronous machine, x,(1) 
P.,(0)<0,P,(0}>0,8(0)>0 
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‘Synchronous machine, x,(t) 
P,(O<0.P,(0)=0,5(0)>0 


Synchronous machine, x(t) 
P,(0}<0.P,(0)=0,5(0)>0 
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Figure 34: 6(0) > 0, 6(0) > 0, P.(0) > 0. 
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Synchronous machine, x,(t) 
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‘Synchronous machine, x,(t) 
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Figure 35: 6(0) > 0, 6(0) > 0, P.(0) > 0. 
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Figure 36: 6(0) > 0, 6(0) > 0, P.(0) > 0. 
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Figure 37: 6(0) > 0, 6(0) > 0, P.(0) > 0. 
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Simulation to study the effect of 5(0) 


The next result was simulated with 


and the following system parameters, 
P,=127, H=0.01, m=21, m=2.5, m=16, 7T=7.7,Erp =1.28, D=0.052. 
Figure 38 shows the result when 


60) ==, 7=0,1,...,6. 


Figure 38 shows that the initial value 6(0) does not have the effect on the steadystate value of 
6(t), further result done on a longer simulation time did show that the steadystate value was 
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‘Synchronous machine, x,(t) 
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Synchronous machine, x,(t) 
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Figure 38: Synchronous Generator Py», = 1.27, P.(0) = 2.3, 6(0) = 


nT 
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Simulation to study the effect of P. 


To study the effect of P., let 
P,=n, n=0,1,2,...,6. 


Figure 39 was simulated with 


together the following system parameters, 
Py = 1.28, H=0.014, m=1.98, m=2.67, 73=1.57, 7=9.1,Erp = 1.37, 


Figure 39 shows that 
T : T 
ri jim 6(é) = 0.75 rad < 5 
Also, 


lim 6 =0 
t-0co 


D = 0.0588. 


except when P.(0) = 0 where 0(t) increased greatly. The steady state values of P.(t) was 


lim P. = 1. 
t>0co 


Further simulations done with the same system and parameters, but with P.(0) < 0 gave unstable 


result, ie 
lim 6(t) > c, 


too 


no matter whether P,,, > 0 or < 0. 
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Figure 39: Synchronous Generator P,, = 1.28, 6(0) =1, P.(0) =n, n=0,1,...,6. 
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Simulation to study the effect of P,, 


Figure 40 shows the result simulated to study the effect of P,,. The system parameters were set to 
H=0.021, ,=2.98, m=2.67, ng=1.74, 7=101,Epp=1.12, D=0.0672, 
while the initial conditions of the states were set to 


6(0)=1, P,(0) = 0.98. 


Figure 40 shows that there is a maximum point max P,, for the system to remain stable, and 
from this simulation 
max P,, & 2.0. 


Po > max 7 
will result in 

lim d6(t) > oo. 

t-00 


The angle 4(t)|, ,., depends on P,,, and d(t)|, ,., increases with P,, until max P,,, is reached, the 
corresponding 6 at this point (d,,) from this simulation was 


max 0, < 0.8 rad. 


Also, P. at steadystate decreased as P,, increased. When P,, was 2.0 and the system became 
unstable lim;_,.. P. decreases to below 0.8. 
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Figure 40: Synchronous Generator Py», = 1.28, 6(0) =1, Pr(0) =n, n = 0,0.4,0.8...,2.0. 
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Figure 41 shows a plot between P,, and the steadystate P., where I assumed that 
kim Pe(b) PE (200). 


Figure 41 shows how P.|, ,.. decreases sharply at the point P,, + 1.9. Also, 


P_v.s. P. 
m e 


for a synchronous generator 
T T T 


a.°0.7+ 4 


0.5 a 


0.47 


Figure 41: P,, v.s. P.(200) of a synchronous generator P.(0) = 0.98, 6(0) = 1, Pr(0) = n, 
= 0.01030. 06.4.0: 


Palicve eo Oa. tor JP 9, 


Technical record #3 by Kittisak Tiyapan, Tokyo Institute of Technology 61 


Figure 42 shows a graph between P,, and 6 of the same synchronous generator as above. 
Figure 42 shows how 6(t) also runs away 


Pvs.5 
m 
‘m VS: for a synchronous generator 
nous generator T T T T 


12000 


10000 - 


8000 - 


© 6000; 


4000 


2000- 


Figure 42: P,, v.s. 6 synchronous generator P.(0) = 0.98, 6(0) = 1, P,(0) =n, n = 
0, 0.03, 0.06,...,4.0. (b) is a closed-up detail of (a) 


The saturation of the graph on the right hand side of Figure 42 could be due to the fact that the 
simulation was done upto t = 200 and not t + oo, according to the approximation mentioned earlier 
in Equation . (Actually this was confirmed by doing another simulation with a longer simulation 
time. When tina) = 400 was used instead of t pinay = 200 the graph of Figure 42 (a) went up to 
6 2.3 x 10* instead of 6 © 1.2 x 104 as in the figure.) 
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Figure 43 shows a graph between P,,, — P. (which is called the accelerating power of the machine) 
and 6 of the same synchronous generator. Figure 43 shows how the accerating power, 


(P_-P,) vs. C) 
for a synchronous generator 
T T T T T T 


Figure 43: (P,—P.) v.s. 6 synchronous generator P.(0) = 0.98, 6(0) = 1, Pr(O) =n, n = 
0, 0.03, 0.06,..., 4.0. 
Po =Py—- P,, (44) 
swings from negative to positive. From this figure, at 
60.7 


a point is reached where the area between the curve of the graph and the P, axis adds up to zero. 
Beyond this point the machine loses synchronism and the speed increases indefinitely. In practice, 
however, this increase speed may be limitted by physical factors. 
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Equal area criterion 


The equal area criterion [AF77] can give the reason why there exists a maximum dma, Which 6(t) 
may not exceed if the system is to remain stable. The equation is derived in the following manner. 
Let P, be called the accelerating (per unit) power, then 


2H d?5 
ae Ph —- Pe = ia 
Wp dt? 

d?6 = WR 

Ges 2H 


Multiply both sides by 2%, 


Integrate both sides, 


2 6 
@ _ Rf pas 
dt 


dé WR 6 
Be eS lee he ds, 4 
Doda 45) 


Equation 45 gives the relative speed with respect to a reference frame (which moves at a constant 
speed.) 
For stability, 


when rotor is not accelerating, a must be zero, and 


when rotor is accelerating, a must oppose rotor motion (relative to the reference frame), and 


the stability conditions is that 


6 
Wrox St. Pa (max) $0, and f/ Padé = 0. (47) 
‘0 
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Discussion 
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Volterra-Lotka ecosystem 
Consider a Volterra-Lotka ecosystem model written in state space form as 


Z, = (a— bx) x (48) 
XQ a (cx, = d) v2, (49) 


where 
Gb Cd ert) pete Re. 


Some knowledge about this system 


The equilibrium point can be found by letting 


z= f(x) =0, 
which leads to 
pit. jesG 
: iy 2 bj 


as the equilibrium point. 
The nature of this equilibrium point can be found by looking at the eigenvalues of the matrix 
A of the system linearized around the equilibrium point itself, that is 


& = Ax + Bu, (50) 


where A is evaluated from the Jacobian matrix at the equilibrium point, ie. 


0 fal 
= eel (51) 
(48) rE : 


a—br. —bx, | 


The eigenvalues of A are then 


which when a,d > 0 leads to the eigenvalues at 
qv ad. 


The eigenvalues are pure imaginaries means that the equilibrium point is a center of trajectories. 
Figure 44 shows the result of simulations done on this system having 
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Phase plane 
T 


x1 and x2 


Time in seconds 


(a) 


Figure 44: A Volterra-Lotka ecosystem model, x2 is shown with a dotted line (a = 3, b=1, c= 2, 
d=1.5), (a) state variables, (b) state plane. 
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With a signum function feedback 


If an input 
u = k sgn (x; + x2) 
was applied to this system in the following manner 
L1 = (a = bx2) Uy (52) 
(53) 


fq = (cx,-—d)xo+u 


the response was that shown in Figure 45 (k = 1), and Figure 46 (k = —1) 
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Figure 45: A Volterra-Lotka ecosystem model with u = sgn (x; + X2), Z2 is shown with dotted line 
2, d= 1.5), (a) shows state variables, (b) state plane 


(A= 3) b= 16 
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(a) 
Figure 46: A Volterra-Lotka ecosystem model with u = —sgn (x; + X2), Lg is shown with dotted line 
1,c=2, d= 1.5), (a) shows state variables, (b) state plane 


(G—350 
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Two signum function inputs 


Incorporating into Equaiton 48 and 49 the constraints (from reality) that 


20) tp 20, (54) 
and also add inputs 
Uy = A, Sgn (O21 rT tg — 01) (55) 
Uy = Kosgn (Vor, + 2 — 02) (56) 
to obtain 
Z1 = (a—bmax (x9, 0)) max (21,0) + Kisgn (0121 + £2) (57) 
Z2 = (cmax (2,0) — d) max (22,0) + Kosgn (Vor + £2). (58) 


Here k1, Ko, V1, V2, 01, and @2 are constants and 
V1, V2, 01, Q2 = 0. 
The following simulation used the set of parameters 
= 24 0 = 14) .6= 73.0 =2 6.07 = 03,07 = 21,01 =— 65, Os =] 98. 
For the feedback gain «; and k2, four values were used, they were the combinations of 
Ky = +5.6;. andks = =48.3. 


The result from the simulations is shown in Figure 47, 48, 49, and 50. The result from simulation 
when only either one of the two &’s was applied was also obtained, and is shown in Figure 51, 52, 
53, and 54. In all cases, the initial conditions were taken as 


(x1(0), z2(0)) = (0.1, 6.9), (2.1, 4.9), (4.1, 2.9), (6.1, 0.9). 


The plots of u,(t) and u2(t) are shown in this order of initial conditions, ie. left > right and then 
top — bottom. 
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Control input u ‘i 
pute Control input u, 
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Figure 47: Two signum inputs Volterra-Lotka system Kk, = 5.6, k2 = 8.3 (a) x(t), (b) x(t), (c) v1 
- (d) max (21,0) (t), (e) max (x2,0) (t), (f) max (21,0) (t) v.s. max (x, 0) (t), (g) ur(t), (h) 
U2 t 
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Figure 48: Two signum inputs Volterra-Lotka system Kk, = 5.6,K9 = —8.3 (a) x(t), (b) x(t), (c) 
ie . (d) max (x1, 0) (t), (e) max (x2, 0) (t), (f) max (x1, 0) (t) v.s. max (22,0) (t), (g) ui(t), 
U2 t 


Technical record #3 by Kittisak Tiyapan, Tokyo Institute of Technology 


maxx 


Na 
08 tl ot 
' ' 
Control input u, (for the four initial conditions used) 
6 6 
4 4 
2 2 
0 0 
2 -2 
-4 -4 
6 ~~ —— = 
rh nn | 05 1 15 
t t 
6 6 
4 4 
2 2 
eH 0 
-2 -2 
-4 -4 
6 6 
0 05 1 15 05 1 15 2 


Figure 49: 
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Two signum inputs Volterra-Lotka system kK, = 
(d) max (2,0) (t), (e) max (x2, 0) (t), (f) max (x, 0) (t) v.s. max (22,0) (t), (g) ur(t), 


—5.6, kg = 8.3 (a) x1(t), (b) xa(t), (c) 
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Figure 50: Two signum inputs Volterra-Lotka system k, = —5.6, kg = —8.3 (a) x1(t), (b) xo(t), (c) 
a - (d) max (x1, 0) (t), (e) max (x2, 0) (t), (f) max (x1, 0) (t) v.s. max (22,0) (t), (g) ur(t), 
U2 t 
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Figure 51: Two signum inputs Volterra-Lotka system Kk, = —5.6,k2 = 0 (a) 11 (t), (b) xo(t), (c) v1 
v.8. L, (d) max (x1,0) (t), (e) max (x2, 0) (t), (f) max (21,0) (t) v.s. max (22,0) (t), (g) ur(t) 
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Figure 52: Two signum inputs Volterra-Lotka system k1 
v.8. L, (d) max (x1,0) (t), (e) max (x2, 0) (t), (f) max (21,0) (t) v.s. max (29,0) (t), (g) ur(t) 


5.6, K2 = 0 (a) x1 (t), (6) xo(t), (c) x1 
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Figure 53: Two signum inputs Volterra-Lotka system Kk, = 0,Kk2 = —8.3 (a) 11 (t), (b) x2(t), (c) x1 
v.8. L2, (d) max (x1,0) (t), (e) max (x2, 0) (t), (f) max (21,0) (t) v.s. max (22,0) (t), (g) u(t) 
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0 1 2 3 4 5 0 1 2 3 4 5 
t t 
10 10 
5 5 
“9 S 0 
-5 -5 
-10 -10 
0 1 2 3 4 5 1 2 3 4 5 
t t 


Figure 54: Two signum inputs Volterra-Lotka system Kk, = 0, k2 = 8.3 (a) x1(t), (b) xo(t), (c) x4 
v.8. £2, (d) max (x1,0) (t), (e) max (x2, 0) (t), (f) max (21,0) (t) v.s. max (22,0) (t), (g) ua(t) 
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Sensitivity analysis 
To do a sensitivity analysis let (see also [Kha96]) 
CHAN. bo) — Bae 
The state variables can be written as 


#(t,A) = 29+ | FS aR hile: 


Take partial derivative with respect to A, 


_ ft lof Of 
x(t, x) — a Fie a(S, r), ATA; X) ze ys x(s, r), X) ds, (59) 
where ax(t, d) a 
_ a(t, Oy _ 
Gy lit) = a and Dn 0 
Differentiate with respect to t, 
apentts A) = A(t, A)xy(t, A) + Bt, A), xa(to, A) = 0, (60) 
where F 
A(t, \)ay(t, 4) = OME) v=x(t,r), B(t,A) = ee, = x(t,). 
Bi 


Let the sensitivity function be S(t) = x(t, Ag) and obtain the sensitivity equation 
From Equation 48 and 49 the Jacobian matrix with respect to the state variables is 


Of _ es —bx1 | 


Ox CXL2 —d+ cx, 
and the Jacobian matrix with respect to the parameters is 


Of... XL, XX 0 0 
Orn | 0 0 Dyey Ly |) 


Let the nominal values of parameters are 


to obtain 


Of 


Ox 


= 1— 22 —X1 | Of 


7 HD) —l + 21 


nominal nominal | 


Let the sensitivity function be 


Ae ae Je B 


Oa Ob Oc Od 


Tq LE LE Xo 


U3 U5 U7 Ug 
s-| 


| | Ox, Ox, Ox, Ox, | 


nominal 
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Then for the purpose of sensitivity analysis, one simulates the system 


v1 = (1 — £2) Uy x1(0) = 10 
Lo = (xy _ 1) v2 x2(0) 20 
£3 = 3 — %2%3 —%X4 + 2X1 x3(0) X30 
Lq = Lol3—Xgt 11 L4 r4(0) 40 
Zp = 5 —@o%5—X%1F——2XX2 %5(0) = Fs 
Le = L2X5 — Xo t+ 11% re(0) = Xoo 
Ly = XL7—Xol7 — 11Lg r7(0) = 2X70 
Zag = U7 —AXgt+X1Fgt+71X2 Fg(0) = Le 
Lg = Ly — LyX — L1L19 r9(0) = Xoo 
Lig = 2g — Lig + X19 — Tq X49(0) = oo 
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Synchronous machine connected to an infinite bus 


79 


When a synchronous machine is connected to an infinite bus! via a transmission line the current 


model is (see also [Tiy98]) 


tp : ; 7 ip 
ie -L (R+wN) 0 i: 
a) = ed ee 
2a tQ 
wy ee ae ae ee ll ae ee, 
é 0 0 0 0 0 1 0 é 
—Ksiny 
piel 0 0 
Kcosy |, 
0 
0 1 0 = 
0 1 —1 
where r, La, and L, in L, R, and N are replaced by 
R=r+Re, Lg=Latle, Ly=Ly+Le, 
R. and L, being resistance and inductance of the transmission line respectively, and 
y=dO-a. 
The flux linkage model is 
Tz=C2z+D, 
where 
T 
r=|[\ Aw AD Ay Aq w 6], 
ag See) eae ate 
1 0 0 0 
0 0 1 
Ee, Le LuaQ LeLuMaq 
a 0 a lq (1 lg ) Igla @) 
0 it 
28) 
0 0 01 


‘An infinite bus is a constant-voltage, constant-frequency bus (source or sink) 


(62) 


(63) 


(65) 


(66) 
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C= 
_R Lup RLup RLIMp L L wLheL 
cea) lal lal —w [1+ (1 He) Eig 0 0 
rplup ate-(4, Lup rplup 0 0 0 0 
lpla lp F lplp 
rpLup rpLup rp ({ — Lup ) 0 0 0 
Tals tale Lp Lp 
R Iu RLM 
Ww [1 Le (1 Lup )) wleLyp wlheLup —F( =; 2) 7m Q 0 0 
lg lg lale lalp q qd q'Q 
0 0 0 rolma ey Gl 0 0 
lglg Ig lo 
—Lmup) _ Lup Iup Img Imuq@ dr _D 
37,12 q 3T; lal q 37; lalp q 37514 q 37; lala qd Tj 
0 0 0 0 0 1 O 
(68) 
and 
| V3Vinfty sin(d — alpha) ] 
UF 
0 
Ds —V3Vin fty cos(é — a) (69) 
0 
Tm 
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Appendix 


Reponses of some 2”/-order systems with simple state equations 
With trigonometric functions in the state equations 


With systems described by an ordinary differential equation 
eb See K., (70) 


where yx is a trigonometric function and & < 0 is a constant. Written in state equation form the 
system is 


Ly = vr (71) 
fo = Fx(21) +k. (72) 


Figure 55 shows the result when x (71) = sinz,, k = 2.7 for Figure 55 (a) and & = —2.7 for 
Figure 55 (b). The result when # — sina = + is similar to that when # + sina = +x, and thus is 
not shown here. Likewise, the result when 7 + cosx = +& is similar to that when % + sinx = +k, 
and therefore in not shown here either. The plots of the state variables start from the initial point 
of (—5, 5). 


x"+¢sin x=k 
wy state plane 


Figure 55: State variables and state planes of (a) &+sinx = 2.7, and (b) &+sinx = —2.7 
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Figure 56 shows the result when y (#1) = tana, « = 2.7 for Figure 56 (a) and & = —2.7 for 
Figure 56 (b). Figure 56 came across singularities while simulating. 


x"stan x=k 
1 state plane x"tan x=-k 
1 state plane 


Oo, — 


Figure 56: State variables and state planes of (a) ¢ + tanx = 2.7, and (b) #+tanz = —2.7 


Figure 57 shows the result when x (#1) = —tan2,, & = 2.7 for Figure 57 (a) and « = —2.7 for 
Figure 57 (b). 


x"-tan x=k 
SEK, state plane 
4 
SVX, 
-5 A \ o—— 
\ 
-6 \ -50 
7 
so -100 
8 ie 
-9 -150 
-10 -200 
M9 2 4 6 8 10 10. ~OO 10 20 30 ~250 
. se Oo 3 4 6 8 10 
t 1 t 
SD state plane SV. % 
5 T ) 7 10 
[ | 10 y 
| -10 — 
5 a ae 
n “ ¥ 
0 | x = eS 
| = 
-30 
| | 0 Ne 
-40 
J / L \ 0 2 4 6 8 10 
5 5 t 
0 2 4 6 8 10 10 0 10 20 30 
a x, 


Figure 57: State variables and state planes of (a) % — tanx = 2.7, and (b) & — tanxz = —2.7 
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Figure 58 shows the result when y (71) = cscx1, k = 2.7 for Figure 58 (a) and & = —2.7 for 
Figure 58 (b). 


x"+c0sec X=-k 


i state plane 
0 
x"4c08e0 x=k 
SV. X, state plane -10 
250 20 
200 15 — 
150 10 nan 
2100 x5 a 
50 Q a 
OL fe 60 
0 2 4 6 19 
& 10 
0 2 4 6 8 10 =10 
: 8.V.% 
10 
SV. x, 
40 20 5 
35 15 
30 ‘ 
10 
25 cad 
ns 1 5 
20 -10 
0 
15 =18: 
10 es 
-20 
: a i 3 7 @ 10 
0 2 4 6 8 10 =10 


Figure 58: State variables and state planes of (a) % + csc x = 2.7, and (b) + cscx = —2.7 


Figure 59 shows the result when y (z,) = —csca1, K = 2.7 for Figure 59 (a) and « = —2.7 for 
Figure 59 (b). 


x"—cosec x=k 


4 state plane x"—cosec x=-k 


state plane 


Figure 59: State variables and state planes of (a) % — cscx = 2.7, and (b) & — cscx = —2.7 
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Figure 60 shows the result when y (71) = secx1, « = 2.7 for Figure 60 (a) and & = —2.7 for 
Figure 60 (b). 


x"s8ec X=-k 
f state plane 


x"sec X=k Or 
4 state plane 


Figure 60: State variables and state planes of (a) % + secx = 2.7, and (b) + secx = —2.7 


Figure 61 shows the result when y (41) = —sec21, Kk = 2.7 for Figure 61 (a) and « = —2.7 for 
Figure 61 (b). 


x"Lsec X=-k 
4 state plane 


x'-sec x=k 
a state plane 


100 
80 
60 
40 
20 
0 
a 
-20 
0 2 
30 
20 
10 
0 
-10 -10 
0 2 4 6 8 10 0 20 40 60 


Figure 61: State variables and state planes of (a) & — secx = 2.7, and (b) % — secx = —2.7 The 
initial condition for the state variable plots in (b) is (3,3) 
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Figure 62 shows the result when y (x71) = cota1, & = 2.7 for Figure 62 (a) and « = —2.7 for 
Figure 62 (b). 


x"+cot X=-k 
paca state plane 


x"+cot x=k 


ae 8.V. x, state plane 4 | \ 
\ I 25 \ 
VA AK ATTA LELT 
al ffl “t\L\I VIVA 
geal, Mey a da paul srt VT EVE VY 
eae ae ac AE VP 
15 | \ / a | | \ | | \ 
i ie A Ee | | 
0 2 4 6 8 10 
eat) 2 4 6 8 10 -20 0 20 40 t 
: 4 SN. state plane 
SN state plane 5 5 
4 20 =a \ )\ 
| — ae 
ANA Aa AI 
“0 a | | | -10 
-2 | | | 0 | \ | -15 
E \ en 
o 2 4 6 8 10 o 10 20 30 40 60 *o 2 4 6 8 10 8 40-80-20. 100 


t x, t x, 


Figure 62: State variables and state planes of (a) + cotx = 2.7, and (b) &% + cotx = —2.7 The 
initial condition for the state variable plots in is (3,3) 


Figure 63 shows the result when y (21) = — cot 21, « = 2.7 for Figure 63 (a) and & = —2.7 for 
Figure 63 (b). (singularities were encountered while simulating Figure 63) 


x"Lcot x=-k 
SeXy, state plane 


x"Lcot x=k 
SV. X, state plane oF 


0 2 4 6 8 10 -10 0 10 20 30 40 50 
t x, 


Figure 63: State variables and state planes of (a) & — cotx = 2.7, and (b) & — cotx = —2.7 The 
initial condition for the state variable plots in is (3,3) 
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Figure 64 shows the result when y (v1) = sin7' 21, « = 2.7 for Figure 64 (a) and k = —2.7 for 
Figure 64 (b). 


Figure 64: State variables and state planes of (a) # +sin~' = 2.7, and (b) ¢+sin7! x = —2.7 The 
initial condition for the state variable plots in is (—0.8,0.8) for (a) and (0.8, 0.8) for (b) 


Figure 65 shows the result when y (x,) = —sin7! 21, & = 2.7 for Figure 65 (a) and « = —2.7 for 
Figure 65 (b). (singularities were encountered while simulating Figure 65) 
For y (1) = cos”! x, there are singularities within the equation and therefore could not simulate. 
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i) state plane 


t x 


Se state plane 


<) 


1 
0 2 4 6 8 10 -8 -6 4 -2 0 2 4 

x 

1 


(a) (b) 


i 

a 
bbe eb 
Bibs the Bieter ok 


Figure 65: State variables and state planes of (a) # — sin~'x = 2.7 topy = 3.99sec, and (b) 


#—-sin7' x = —2.7 topy = 4.81 sec The initial condition for the state variable plots in is (—0.3, —0.3) 


for (a) (3,3) for (b) 
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Figure 66 shows the result when y (z,) = tan”! 21, Kk = 2.7 for Figure 66 (a) and « = —2.7 for 
Figure 66 (b). 


o & k & b&b 4 o 
x, 

i i fo ae 32 

& > kw oN & Oo @ 


1 
Pen ape sia eo 
bo kW ON AO @ 


° 
bo kh oD Ao @ 
Xp 

° 
bo kW ON AO @ 
=O 


Figure 66: State variables and state planes of (a) %+tan~'x = 2.7 tepy = 0.64sec, and (b) 
¢+tan7!x = —-2.7 tepy =0.61sec The initial condition for the state variable plots in is (—3, —3) 


for (a) and (3, —2) for (b) 

Figure 67 shows the result when y (%,) = —tan7!2,, « = 2.7 for Figure 67 (a) and k = —2.7 
for Figure 67 (b). (integration tolerances were violated while simulating Figure 67) 
With polynomial of state variable terms in the state equations 


Consider Equation 71 and 72 again but this time let 
x(a) = 2", n=—5,-4,,-3,...,5. 


Figure 68 shows the result when y (%1) = 71, « = 2.7 for Figure 68 (a) and k = —2.7 for Figure 
68 (b). 

Figure 69 shows the result when y (x1) = —a1, « = 2.7 for Figure 69 (a) and « = —2.7 for 
Figure 69 (b). (integration tolerances were violated while simulating Figure 69) 
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a” 1 
x"-tan”! x=k x"-tar”! x=-k 


a5 G state plane = SVX, state plane 
100 0} —___— 
x 50 x -50 
o-———_— -100 
-50 -150 
2 4 6 8 10 0 2 4 6 8 10 
t qt: 
mane state plane SM% state plane 
25 10 5 
20 18 be 5 
15 i) 
5 0 
10 -5 
aa aa a ~~ 
5 ; -10 5 
0 Mg -15 
a 
-5 -20 -10 = 
-5 
-10 -25 
2 4 6 8 10 -5 0 5 10 2 4 6 8 10 -10 -5 0 i} 
t x, t x 


Figure 67: State variables and state planes of (a) % — tan~+x = 2.7 topy = 0.75sec, and (b) 
% —tan7!x = —2.7 topy = 0.88sec The initial condition for the state variable plots in is (4, —4) 


for (a) (—4,4) for (b) 


S.V. X state plane ox eo 10 
8000 2 
-0.5 
6000 
-1 
4000 ‘2 
e A156 
2000 -2 
6 ees -2.5 
-2000 7G 2 4 6 8 10 
2 4 6 8 10 
‘ t 
S.v. 
SV % state plane x10" Nt 
8000 10 oP 


oo 
4000 3 4 a 
x “0 it x 0 
15 
2000 © 
-2 —o 
; a 5 -5 —o 
-25 
-2000 10 -3 ~10 
0 2 4 6 8 10 =10 5 0 5 10 0 2 4 6 8 10 10 5 0 5 10 
t x, t x 


(a) (b) 


Figure 68: State variables and state planes of (a) +x = 2.7 topy =0.41sec, and (b) +2 = —2.7 
tcpy = 0.38sec The initial condition for the state variable plots in is (—1,—1) for (a) and (—1,1) 


for (6) 
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hye X"-x=-k 
S.V. X, aie plane SMX, state plane 
8 10 0 10 
-1 
6 7 5 
4 -2 
$f NO om 0 
2 -3 
0 = -4 a 
- & 5 10 
25 2 4 6 8 10 10 0 2 4 6 8 10 
t t 
SVX S.V. X 
4 10 2 10 
2 5 1 5 
“0 “0 x0 “0 
2 5 1 -5 
4 -10 2 10 
0 2 4 6 8 10 5 0 5 10 0 2 4 6 8 10 


Figure 69: State variables and state planes of (a) #-—x = 2.7 tepy = 0.41 sec, and (b) &-x = —2.7 
topu = 0.38sec The initial condition for the state variable plots in is (—1,1) for (a) (—1,-1) for 


() 


Technical record #3 by Kittisak Tiyapan, Tokyo Institute of Technology 91 


Figure 70 shows the result when y (x1) = x7, « = 2.7 for Figure 70 (a) and « = —2.7 for Figure 
70 (b). When # + x? = 2.7 there are singularities within the solution and could not simulate. The 
same thing also happened with 7 — x? = —2.7, therefore neither of the two cases could be simulated. 
When # + 2? = —2.7 simulation (Figure 70 (a)) showed that the response is stable approximately 
in the range of 
—3.28 < x, < 1.61, 


outside of which there appeared singularities in the solution. Likewise Figure 70 (b) and additional 
simulations showed the range of stability to be 


—1.64 < x; < 3.28 


for % — x? = 2.7. 


x" ok x" =k 
1 state plane prise: 3 state plane 


Figure 70: State variables and state planes of (a) ¢+x* = —2.7 tcpy = 0.13sec, and (b) ¢-x? = 2.7 
tcpu = 0.14sec 
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When # + 2° = +2.7 there are singularities within the solution and could not be simulated. 
Figure 71 shows the result when % — 73? = +2.7 


x'aak x’ =—k 


4 state plane SV. X; state plane 


Mal state plane SS 


Figure 71: State variables and state planes of (a) i—x° = 2.7, and (b) #2? = —2.7 tepy = 0.15 sec 
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